使用半方差模型

# 如何理解半方差图：
https://pro.arcgis.com/zh-cn/pro-app/help/analysis/geostatistical-analyst/understanding-a-semivariogram-the-range-sill-and-nugget.htm

## 如何利用渔网稀疏解决空间自相关的问题：
https://kevintshoemaker.github.io/NRES-746/SDM_v6.html#species:_common_chuckwalla_(sauromalus_ater)

## 如何在宏观尺度上计算半变异函数：
https://www.aspexit.com/en/implementing-variograms-in-r/
https://www.r-spatial.org/r/2016/02/14/gstat-variogram-fitting.html
# 参见：
library(sp)
data(meuse)
head(meuse)
# no trend:
coordinates(meuse) = ~x+y
variogram(log(zinc)~1, meuse)
# residual variogram w.r.t. a linear trend:
variogram(log(zinc)~x+y, meuse)
# directional variogram:
aa <- variogram(log(zinc)~x+y, meuse, alpha=c(0,45,90,135))
bb <- variogram(log(zinc)~1, meuse, width=90, cutoff=1300)
plot(bb)
variogram(Chuckwalla ~1, k.chuck.trn, cutoff = 50000, xlab= "Distance (m)", ylab="Semivariance")

## 半方差图中，半方差与空间距离之间存在空间依赖关系；
# 一般半方差接近于平行分布时表示无空间自相关；
# 而一般半方差曲线会在较短距离时倾斜分布，但存在一个问题。如果存在一个较强的空间自相关，则这种趋势会延续到更远的距离，一般以米为单位；当存在较强的空间自相关关系时，解决的办法是采用渔网稀疏方法，设定单元格密度范围内进行空间稀疏，保留指定数量的分布点。再进行空间拟合，这时候会看到半方差的倾斜值对应的x轴会降低，从而仅保留较低尺度上的 空间依赖关系，从而减弱了空间自相关的存在度。

## 在sdm中解决空间自相关的方法：
## 同样基于半方差的方法：
# 参考：http://rstudio-pubs-static.s3.amazonaws.com/9687_cc323b60e5d542449563ff1142163f05.html

第一种就是空间稀疏吗，但即使是空间稀疏可能检测到残余空间自相关(RSC、ASC)。空间稀疏可能会丢失物种分布信息；
第二种方式：使用统计方法，这些方法允许我们考虑（甚至消除）空间自相关。我将要介绍的所有方法都假设空间平稳，这意味着空间自相关和其他协变量的影响在整个区域中都是恒定的。这种假设有时可能会出现问题。Dormann等。提供了以下示例：如果空间自相关的主要原因是分散的，则如果动物移到或多或少受到限制的不同栖息地，就会破坏平稳性。不幸的是，很少有方法可以解决非平稳性问题。

# 第二种方法具体来讲：
这种情况下的最佳解决方案是找到缺失的协变量并将其包含在模型中。不幸的是，这并非总是可能的，我们必须找到另一种解决方案。一种可能性是对该剩余空间自相关建模，并将此模型包括在我们的线性丰度模型中。

## 工作方法：
为了说明空间自相关，我们将使用软件包nlme中提供的名为GLS（广义最小二乘）的过程。实际上，GLS使用参数函数直接在方差-协方差矩阵中对空间协方差结构进行建模。但是首先我们将忽略空间自相关，并使用gls函数（而不是lm）重新拟合介绍中的模型。结果将是相同的，但是稍后使用AIC进行模型比较时，我们将需要此模型（即，我们无法使用lm将模型拟合的AIC与使用gls进行拟合的AIC进行比较）。
## 代码如下：
# 先构建一个空间回归模型，并计算半方差：
在方差图上，我们看到半方差明显随距离而增加。我们确认残差中存在空间自相关。
m2 <- gls(log1p(counts) ~ elevation + I(elevation^2), data = sampledata)
vario2 <- Variogram(m2, form = ~x + y, resType = "pearson")
plot(vario2, smooth = TRUE, ylim = c(0, 1.2))

## 检测到空间残差中存在空间自相关之后，我们使用不同的相关结构拟合模型，然后使用AIC选择最佳模型（我们还将比较没有相关性的初始模型）。
m3 <- gls(log1p(counts) ~ elevation + I(elevation^2), correlation = corExp(form = ~x + 
    y, nugget = TRUE), data = sampledata)
m4 <- gls(log1p(counts) ~ elevation + I(elevation^2), correlation = corGaus(form = ~x + 
    y, nugget = TRUE), data = sampledata)
m5 <- gls(log1p(counts) ~ elevation + I(elevation^2), correlation = corSpher(form = ~x + 
    y, nugget = TRUE), data = sampledata)
m6 <- gls(log1p(counts) ~ elevation + I(elevation^2), correlation = corLin(form = ~x + 
    y, nugget = TRUE), data = sampledata)
    m7 <- gls(log1p(counts) ~ elevation + I(elevation^2), correlation = corRatio(form = ~x + 
    y, nugget = TRUE), data = sampledata)

AIC(m2, m3, m4, m5, m6, m7) # 这里应该是m3 的AIC值最低；
绘制拟合的方差图以检查拟合是否正确。
vario3 <- Variogram(m3, form = ~x + y, resType = "pearson")
plot(vario3, smooth = FALSE, ylim = c(0, 1.2))
作为最后的检查，我们绘制归一化残差的样本变异函数图。这些残差是标准化残差预乘以估计的误差相关矩阵的平方根平方的倒数。因此，如果适当考虑了剩余空间自相关，我们就不会在新的变异函数图中看到任何趋势。
vario4 <- Variogram(m3, form = ~x + y, resType = "normalized", maxDist = 30)
plot(vario4, smooth = FALSE, ylim = c(0, 1.2))

## 这里是基于预测结果重新数据：
resp <- expm1(predict(m3, newdata = data.frame(elevation = values(elev), x = coordinates(elev)[, 
    1], y = coordinates(elev)[2])))
resp <- rasterFromXYZ(cbind(coordinates(elev), resp))
plot(resp, main = "Predicted abundances")

## 关于解决空间残差的方法：（也就是模型未能拟合部分）：作者提出了一些思路：
通用克里金法在数学上等效于以下两步过程：首先使用协变量将线性模型拟合到您的数据，然后对残差执行普通克里金法，然后将克里金残差添加到线性模型的预测中。您将获得的新预测与通用克里金过程提供的预测相同。这种计算通用克里金法预测的方法通常称为回归克里金法。我们立即看到了这种替代参数化的主要优点之一：可以使用例如泊松分布将GLM（广义线性模型）拟合到我们的计数数据中；然后，我们可以使用普通克里金插值GLM残差并将其添加到我们的预测中。当然，我们应该首先在链接规模上将kriged残差添加到预测中，然后将结果反向转换。不幸的是，使用这种方法获得相应的标准误差要复杂得多。

### 计算莫兰指数：
### 莫兰指数进行运算 ####
install.packages("lctools")
library(lctools)
data(GR.Municipalities)
Coords<-cbind(GR.Municipalities@data$X, GR.Municipalities@data$Y)
#using an adaptive kernel
bws <- c(3, 4, 6, 9, 12, 18, 24)
moransI.v(Coords, bws, GR.Municipalities@data$Income01)